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Abstract 

Store-induced  limit-cycle  oscillation  of  a rectangular  wing  with  tip  store  in  transonic  flow  is  simulated 
using  a variety  of  mathematical  models  for  the  flow  field:  transonic  small-disturbance  theory  (with  and 
without  inclusion  of  store  aerodynamics)  and  transonic  small-disturbance  theory  with  interactive 
boundary  layer  (without  inclusion  of  store  aerodynamics).  For  the  conditions  investigated,  limit-cycle 
oscillations  are  observed  to  occur  as  a result  of  a subcritical  Flopf  bifurcation,  and  are  obtained  at  speeds 
lower  than  those  predicted  ( 1 ) nonlinearly  for  clean- wing  flutter,  and  (2)  linearly  for  wing/store  flutter. 
The  ability  of  transonic  small-disturbance  theory  to  predict  the  occurrence  and  strength  of  this  type  of 
limit-cycle  oscillation  is  compared  for  the  different  models.  Solutions  computed  for  the  clean  rectangular 
wing  are  compared  to  those  computed  with  the  Euler  equations  for  a case  of  static  aeroelastic  behavior 
and  for  a case  of  forced,  rigid-wing  oscillation  at  Mach  0.92. 

Nomenclature 

a angle  of  attack,  degrees 

c chord,  ft 

eg,  ea  center  of  gravity  and  elastic  axis,  ft  from  leading  edge 
CL,,  Cm  lift  and  moment  coefficients  (moment  about  leading  edge) 
i,j,  k computational  index  coordinates 

/ span,  ft 

M Mach  number 

Re  Reynolds  number  (based  on  root  chord) 

p density,  slugs/ft’ 

t time,  nondimensional  (based  on  freestream  velocity  and  wing  chord) 

T temperature,  °R 

X thickness  to  chord  ratio 

U velocity,  ft/sec 

x,  y,  z physical  coordinates  (streamwise,  spanwise,  vertical),  feet 

C,  structural  damping  coefficient,  nondimensional 

Subscripts 

°°  freestream  condition 

.v,  w store  or  wing  property,  respectively 

o wing  root  property 


Paper  presented  at  the  RTOAVT  Symposium  on  " Reduction  of  Military  Vehicle 
Acquisition  Time  and  Cost  through  Advanced  Modelling  and  Virtual  Simulation  ", 
held  in  Paris,  France,  22-25  April  2002,  and  published  in  RTO-MP-OS9. 
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Introduction 

High-performance  fighter  aircraft  with  external  stores  are  required  to  operate  with  high  manoeuvrability 
in  the  transonic  flight  regime.  In  this  regime,  the  potential  exists  for  encountering  transonic  nonlinear 
flutter,  known  as  limit-cycle  oscillation  (LCO).  LCO  is  a limited  amplitude,  self-sustaining  oscillation 
produced  by  an  aerodynamic-structural  interaction,  which  for  the  cases  of  interest,  is  exasperated  by  the 
occurrence  of  shock  waves  on  the  surface  of  the  wing  and/or  stores.  LCO  results  in  an  undesirable 
airframe  vibration  and  limits  the  performance  of  the  flight  vehicle. 

The  main  goal  of  the  current  work  is  to  determine  the  range  of  applicability  of  models  of  varying  fidelity 
to  the  numerical  prediction  of  store-induced  LCO.  This  form  of  LCO  typically  occurs  near  linear  flutter 
boundaries  in  the  nonlinear,  transonic  regime  (Mach  number  ranging  between  0.8  and  1.1),  suggesting 
that  classical  flutter  predictions  using  linear  aerodynamic  theories  can  be  applied  to  the  identification  of 
lightly  damped  modes  that  may  nonlinearly  participate  in  LCO.  Indeed,  using  traditional  approaches, 
Denegri  (2000)  had  limited  success  in  relating  observed  store-induced  LCO  to  “hump”  (or  “soft” 
crossing)  modes  found  in  velocity-damping  diagrams.  However,  in  many  cases,  the  linear  approach  is 
inadequate  in  predicting  response  characteristics  of  vehicle  configurations  in  the  transonic  regime. 

The  transonic  regime  differs  from  the  linear,  subsonic  regime  by  the  appearance  of  shocks.  These 
structures  may  strongly  interact  with  vehicle  boundary  layers,  with  the  possible  consequences  of  flow 
separation  or  significant  shock  movement.  Following  advances  in  nonlinear  modeling  and  computer 
hardware,  nonlinear  aeroelastic  predictions,  including  viscous  effects  and  manoeuvre  loads,  are  tractable 
for  reasonably  complex  configurations  (Melville  (2001,  2002),  Farhat  el  a!.  (2000)).  Still,  the  use  of 
modeling  techniques  that  account  for  viscosity  is  too  computationally  demanding  for  preliminary  design. 
In  a coordinated  manner,  we  examine  the  ability  of  aeroelastic  models  of  varying  fidelity  to  predict 
accurately  LCO  onset  and  amplitude.  Models  based  on  linear  analysis,  transonic  small-disturbance 
theory  (TSDT),  and  TSDT  with  interactive  boundary  layer  are  considered.  Through  this  approach,  we 
discern:  (1)  the  limitations  of  linear  theory  for  LCO  prediction  vis-a-vis  the  simplest  nonlinear  theory 
capable  of  producing  weak  shocks;  (2)  the  ability  of  TSDT  to  predict  store-induced  LCO  in  inviscid  flow, 
and  (3)  the  effects  of  viscosity  on  store-induced  LCO.  This  work  provides  the  increased  understanding  of 
the  LCO  phenomenon,  while  also  serving  to  determine  the  range  of  applicability  and  computational  cost 
of  various  modeling  techniques  to  the  prediction  of  this  phenomenon. 

Three  computational  methodologies  are  employed  in  this  investigation:  the  MSC/NASTRAN  aeroelastic 
analysis  program,  the  TSDT-based  NASA/LaRC  CAPTSDv  computational  aeroelasticity  algorithm  for 
inviscid  and  viscous  flow,  and  the  AFRL  ENS3DAE  Euler/Navier-Stokes  computational  aeroelasticity 
algorithm.  MSC/NASTRAN  is  used  in  the  development  and  analysis  of  structural  models,  the  prediction 
of  linear  aeroelastic  response  of  selected  wing/store  configurations,  and  the  search  for  configurations 
speculated  herein  to  produce  LCO.  The  CAPTSDv  algorithm  is  used  to  carry  out  relatively  fast 
aeroelastic  analyses  for  inviscid  flow  and  viscous  flow  through  interactive  boundary  layer  coupling, 
starting  with  cases  previously  identified  by  MSC/NASTRAN  as  LCO-susceptible.  ENS3DAE  is  then 
applied  to  validate  CAPTSDv  predictions  for  static  aeroelastic  and  dynamic  rigid  behaviour  for  a 
representative  set  of  flow  conditions. 

Problem  Formulation 

The  wing  studied  herein  is  derived  from  the  “heavy”  version  of  the  original  Goland  wing.  Like  the 
original,  the  heavy  wing  is  structurally  represented  by  a beam,  but  with  additional  non-structural  mass, 
as  defined  by  Eastep  and  Olsen  (1980).  This  latest  version,  referred  to  as  the  Goland1  wing,  is  a heavy 
wing  modeled  with  a box  structure  to  enable  a variety  of  store  attachment  options. 

Geometry 

The  Goland'  wing  is  rectangular  and  cantilevered  from  an  infinite  midplane.  A planform  schematic  is 
given  in  Figure  la  (including  store)  and  geometric  parameters  are  assigned  values  in  Table  1.  The  airfoil 
section  is  assumed  to  be  constant  over  the  spanwise  extent  of  the  wing  and  is  chosen  to  be  that  of  a 
symmetric,  parabolic-arc  airfoil,  defined  by  z = 2 Tw  x ( 1 - x/cw)  (0  < x < c,„).  The  wing-tip  store  is 
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mounted  flush  to  the  wing  tip  (see  Figure  la).  The  sectional  shape  of  the  store  is  also  assumed  to  be 
described  by  a parabolic  arc  uniform  in  the  spanwise  direction.  The  leading  edge  of  the  store  is  offset 
(defined  positive)  upstream  of  the  wing  leading  edge,  thereby  providing  the  shape  formula: 

Z = 2 Tx  (x  + Coff)  [1  - (x  + Coff)/cs]  (- Caff  <X<CS-  Coff). 

Wing  Structure  and  Inertial  Properties 

The  wingbox  structure  finite  element  model  is  built-up  from  shear  panels,  modelling  the  spars  and  ribs, 
and  membrane  elements,  modelling  the  wing  skins.  In  addition,  rod  elements  are  included  to  model  spar 
and  rib  caps  as  well  as  posts  that  connect  the  wing  skins  at  every  spar/rib  intersection.  Every  element  is 
modelled  using  the  same  material  properties  (shown  in  Table  1).  Notice  that  the  density  of  this  “material” 
is  extremely  low.  The  model  was  built  this  way  to  allow  the  stiffness  and  mass  properties  to  be 
decoupled.  This  decoupling  allows  the  elements  to  be  sized  to  tune  the  model  to  match  the  structural 
dynamic  characteristics  of  the  beam  model  of  the  heavy  Goland  wing. 


Parameter 

Value 

c l 

6,  20 

Cs,  ls 

10,  1 

coff 

3 

ea 

2 

tw,  Tv 

0.04,  0.036 

Young’s  Modulus 

1.4976  x 109  slugs/ft2 

Shear  Modulus 

5.616  x 10s  slugs/ ft2 

Structural  Density 

0.0001  slugs/ft3 

Table  1:  Selected  values  of  geometric  and  material  parameters. 


The  geometry  of  the  wing  is  simple.  The  origin  is  at  the  mid-height  of  the  root  of  the  leading  edge  spar. 
The  three  spars  are  un-swept  and  placed  at  0,  2 and  4 ft  along  the  positive  x-coordinate.  The  eleven  ribs 
are  evenly  spaced  on  2 ft  centres  along  the  positive  y-coordinate.  The  shear  elements  are  defined  by  the 
intersections  of  the  spars  and  ribs,  with  10  elements  per  spar  and  2 elements  per  rib.  Next,  each  spar  and 
rib  is  0.33334  ft  high  and  each  cell  defined  by  the  spars  and  ribs  is  capped  with  a single  wing  skin 
membrane  element.  This  results  in  a total  of  40  skin  elements.  Finally,  rods  are  added  on  the  top  and 
bottom  of  each  shear  element  and  at  every  spar/rib  intersection.  The  total  number  of  rod  elements  is  137, 
consisting  of  60  spar  caps,  44  rib  caps  and  33  posts. 

The  mass  properties  of  this  wing  are  modelled  by  placing  lumped  masses  with  no  rotational  inertia  at 
each  grid  point.  The  lumped  masses  are  sized  to  match  the  mass  properties  (total  mass,  eg,  and  inertia)  of 
the  heavy  Goland  wing.  For  the  beam  model  Goland  wing,  the  mass  properties  are  simply  modelled  as 
lumped  mass  with  rotational  inertia  centred  at  the  eg  of  sections  placed  at  2 ft  intervals.  For  the  internal 
rib  locations  of  the  (i.e.,  2 ft,  4 ft,  6 ft,  etc.)  beam  model,  the  lumped  mass  properties  were  a mass  of 
22.496  slugs,  with  a rotational  inertia  in  pitch  of  50.3396  slug-ft2,  centred  at  x = 2.6  ft.  For  the  external 
rib  locations  (i.e.,  0 ft  and  20  ft),  the  lumped  mass  properties  were  halved.  In  the  built-up  wing  model,  the 
masses  used  for  the  internal  rib  locations  are  as  follows:  1.9650  slugs  at  each  leading  edge  point,  3.9442 
slugs  at  each  point  on  the  centre  spar;  and  5.3398  slugs  for  the  trailing  edge  points.  Also,  similar  to  the 
beam  model,  the  masses  used  at  the  external  spar  locations  are  half  of  the  amount  used  at  the  internal 
locations. 
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The  final  step  in  developing  the  built-up  model  is  sizing  the  elements  so  that  its  structural  dynamic 
characteristics  match  those  of  the  beam  model  Goland  wing.  For  this  model  development,  the  elements 
are  sized  to  minimize  the  error  between  the  first  three  frequencies  of  the  built-up  and  beam  models.  The 
boundary  condition  for  each  model  is  cantilevered.  To  maintain  symmetry,  the  elements  are  grouped  by 
their  component.  The  thicknesses  (in  ft)  for  the  2-D  elements  are:  upper  and  lower  wing  skins  - 0.0155; 
leading  and  trailing  edge  spar  - 0.0006;  centre  spar  - 0.0889,  and  rib  - 0.0347.  For  the  1-D  elements,  the 
areas  (in  ft2)  are:  posts  - 0.0008;  leading  and  trailing  edge  spar  caps  - 0.0416;  centre  spar  cap  - 0.1496, 
and  rib  cap  - 0.0422. 

Store  Mass  and  Linkage 

The  store  configuration  examined  in  this  work  is  that  of  a tip  store.  This  store  structure  is  modelled  as  a 
series  of  rigid  bar  elements  that  result  in  a 10  ft  long  rigid  bar.  The  resultant  bar  is  centred  0.5  ft  outboard 
of  the  wing  tip  and  2 ft  aft  of  the  wing  leading  edge.  The  store  is  then  rigidly  connected  to  the  six  wing 
tip  grid  points.  This  is  accomplished  by  defining  a seventh  wing  tip  grid  point  at  the  centre  of  the  rib 
structural  model  (i.e.,  x = 2 ft,  y = 20  ft,  z = 0 ft)  and  connecting  this  point  to  the  other  six  points  with  a 
MSC/NASTRAN  RBE3  element.  This  results  in  the  displacement  of  the  seventh  grid  point  being  the 
average  of  the  displacements  of  the  original  six  points.  Finally,  an  RBAR  element  is  used  to  connect  the 
seventh  wing  tip  point  to  the  corresponding  point  (i.e.,  x = 2 ft,  y = 20.5  ft,  z = 0 ft)  of  the  tip  store  model. 

The  mass  properties  of  the  tip  store  are  chosen  to  match  the  properties  of  one  section  of  the  wing:  a mass 
of  22.496  slugs  and  a rotational  inertial  of  50.3396  slug-ft2.  During  this  study,  the  position  of  the  store 
mass  is  fixed  at  the  lateral  and  vertical  centres  of  the  store  (i.e.,  y = 20.5  ft  and  z = 0 ft)  and  varied  in  the 
streamwise  direction. 

Computational  Aeroelasticity  Methods 

CAPTSD  and  CAPTSDv 

CAPTSD  solves  the  three-dimensional,  transonic,  small-disturbance,  potential-tlow  equations  for  partial 
and  complete  aircraft  configurations  (Batina  (1988,1989)).  The  standard,  most  widely  distributed  version 
of  the  program  computes  inviscid,  compressible  flow  for  combinations  of  wings,  fuselage,  horizontal  tail, 
bodies/stores,  and  rectangular  planform  vertical  surfaces.  CAPTSD  solves  the  aerodynamic  equations  of 
motion  using  a time-accurate  algorithm  that  is  capable  of  simulating  both  steady  and  unsteady  flow 
(Batina  (1992)).  The  method  is  capable  of  computing  aeroelastic  interactions  by  coupling  the 
aerodynamic  module  with  a structural  dynamics  simulation.  The  structural  dynamics  of  horizontal 
surfaces  are  simulated  in  CAPTSD  by  using  a modal  structural  model.  The  structural  analysis  is  coupled 
to  the  aerodynamic  analysis  by  a process  that  transfers  generalized  aerodynamic  forces  and  generalized 
displacements  between  the  aerodynamics  and  structural  dynamics  modules.  Using  this  approach, 
CAPTSD  is  capable  of  simulating  both  static  and  dynamic  aeroelastic  phenomena. 

A viscous-inviscid  interaction  version  of  CAPTSD,  known  as  CAPTSDv,  has  been  developed  (Flowlett 
(1987),  Edwards  (1993))  and  applied  to  a variety  of  problems  involving  mildly  separated  and  separation 
onset  flows  (Edwards  (1998)).  The  method  couples  the  inviscid  CAPTSD  algorithm  with  an  inverse 
integral  boundary  layer  model.  The  boundary  layer  equations  are  solved  in  a quasi-steady  formulation 
similar  to  that  recommended  by  Green  et  al.  (1977).  The  outer  inviscid  solution  and  the  inner  viscous 
solution  are  computed  independently  and  are  coupled  using  an  active  control  mechanism  that  minimizes 
coupling  errors  for  unsteady  flows.  In  this  investigation,  a single  version  of  the  CAPTSDv  algorithm  is 
used  to  perform  both  inviscid  and  viscous  aeroelastic  analysis  (execution  mode  controlled  by  user). 

CAPTSD  solves  the  equations  of  motion  on  a sheared  Cartesian  grid  system  where  lifting  surfaces  are 
modeled  as  thin  plates.  Thickness  and  camber  information  for  the  upper  and  lower  surfaces  of  each 
lifting  surface  is  supplied  through  a set  of  surface  slopes  that  are  specified  as  boundary  conditions  for  the 
algorithm.  Similarly,  structural  mode  shapes  are  supplied  as  surface  slope  perturbations.  This  approach 
greatly  simplifies  the  modelling  task  required  for  an  aeroelastic  analysis,  since  the  grids  are  typically 
simple  to  generate,  and  no  moving  grid  algorithm  is  required  for  the  aeroelastic  simulation.  CAPTSD  is 
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1-1.5  orders  of  magnitude  more  computationally  efficient  than  higher  fidelity  methods,  and  offers  the 
potential  for  nonlinear  aerodynamics  analysis  within  a design  framework  using  the  CAPTSD 
methodology. 

ENS3DAE 

Aeroelastic  analysis  of  the  Euler  and  Navier-Stokes  equations  is  carried  out  with  the  Euler/Navier-Stokes 
3-Dimensional  Aeroelastic  (ENS3DAE)  method.  The  Lockheed-Georgia  Company,  under  contract  to  the 
Air  Force  Wright  Laboratory,  developed  ENS3DAE  in  the  late  1 980's  (Schuster  et  al.  (1990)).  This 
program  has  been  used  to  solve  numerous  aerodynamic  and  aeroelastic  problems  about  a wide  range  of 
geometries  including  wings,  wing/fuselage,  wing/control  system,  propulsion,  and  integrated 
airframe/propulsion  configurations  (Smith  et  al.  ( 1 996),  Schuster  et  al.  (1998),  Lewis  and  Smith  (2000)). 

ENS3DAE  solves  the  full  three-dimensional  compressible  Reynolds  averaged  Navier-Stokes  equations, 
the  thin  layer  approximation  to  these  equations,  or  the  Euler  equations  using  an  implicit  approximate 
factorization  algorithm.  Central  finite  differences  are  used  for  spatial  discretization,  and  a three- 
dimensional  implementation  of  the  Beam-Warming  implicit  scheme  is  employed  for  the  temporal 
integration.  Blended  second-  and  fourth-order  dissipation  is  added  to  the  explicit  right-hand-side  of  the 
equations,  and  implicit  second-order  dissipation  is  added  to  improve  the  diagonal  dominance  of  the 
matrix  system.  For  time-accurate  cases,  global  or  local  time  stepping  within  an  inner  sub-iteration  loop 
for  each  physical  time  step  is  employed.  The  sub-iteration  procedure  effectively  removes  the  stability 
limit  on  the  time-step  for  unsteady  flow  cases,  and  improves  the  temporal  accuracy  of  the  method, 
provided  that  relevant  time-scales  within  the  computed  flow  field  are  captured  by  the  selected  time  step 
and  that  a sufficient  number  of  sub-iterates  are  computed. 

The  method  accepts  either  single-  or  multiple-block  curvilinear  grids.  Boundary  conditions  are  imposed 
explicitly  on  each  computational  face  of  each  grid  block  and  the  current  release  of  the  program  requires  a 
one-to-one  match  of  grid  points  at  block  interfaces.  Turbulence  characteristics  are  predicted  using  the 
Baldwin-Lomax  algebraic  turbulence  model  or  the  Johnson-King  model  (Huttsell  et  al.  (2001)).  The 
code  is  written  to  take  advantage  of  vectorization;  directives  for  parallel  operation  on  shared  memory 
processors  are  also  included  in  the  programming.  The  method  is  regularly  executed  on  eight  or  more 
processors. 

A linear  generalized  mode  shape  structural  model  is  closely  coupled  with  the  aerodynamic  method  to 
analyse  structurally  flexible  vehicles.  ENS3DAE  uses  a highly  efficient,  grid-motion  algorithm  for 
aeroelastic  and  control-surface  simulations  that  is  based  on  an  algebraic  shearing  technique.  Since 
dynamic  aeroelastic  and  oscillating  control  surface  simulations  require  grid  models  that  deform  in  time, 
the  algorithm  now  enforces  the  Geometric  Conservation  Law  (Thomas  and  Lombard  (1979)). 

Results 

Physical  Conditions 

The  aeroelastic  analysis  is  carried  out  with  enforced  consistency  between  velocity  and  dynamic  pressure, 
assuming  constant  density  at  sea-level  conditions.  Mach  number  and  Reynolds  number  are  treated  as 
independent  variables,  such  that  match-pointed  conditions  are  not  achieved.  Reynolds  number  is  not 
varied  in  this  investigation.  Structural  damping  is  assumed  to  vanish  for  all  baseline  cases  investigated. 
The  selected  values  of  various  parameters  are  summarized  in  Table  2. 

Summary  of  Grid  Construction 

Three  grids  are  constructed  for  the  CAPTSDv  calculations  reported  in  this  paper.  Owing  to  the  geometry 
of  the  wing/store  configuration  and  the  mid-plane  formulation  of  the  surface  boundary  condition  used  in 
CAPTSDv,  these  grids  are  rectilinear.  Clustering  of  grid  points  is  enforced  along  the  edges  of  the 
geometry  and  normal  to  the  wing  and  store  surfaces.  Grids  are  generated  by  first  computing  grid-point 
distributions  in  each  of  the  three  coordinate  directions  external  to  CAPTSDv,  subject  to  the  specified 
clustering  conditions,  and  then  using  the  generator  internal  to  CAPTSDv  to  obtain  the  full  rectilinear 
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grid.  The  first  grid  (Gl)  is  used  for  inviscid,  clean-wing  and  wing/store  (store  mass  only)  computations; 
the  second  grid  (G2)  is  used  for  inviscid,  wing/store  computations,  and  the  third  grid  (G3)  is  used  for 
viscous  clean-wing  and  wing/store  (store  mass  only)  computations.  The  values  of  parameters  governing 
grid  construction  are  given  in  Table  3. 


Parameter 

Value 

A- 

0.0023771 

7L 

518.67 

Re^ 

15  x 106 

a„ 

0 

<r 

0 

Table  2:  Selected  values  of  physical  parameters. 


Parameter 

Grid  Gl 

Grid  G2 

Grid  G3 

Domain  Size  ( x,y ,z) 

(41,13.5,40) 

(41,13.5,40) 

(41,13.5,40) 

Grid  Dimensions  (ij, k) 

(90,55,60) 

(90,55,60) 

(150,55,60) 

Leading-Edge  Spacing 

0.01 

0.01 

0.008 

Trailing-Edge  Spacing 

0.01 

0.01 

0.008 

Normal  Spacing  (top,  bottom) 

(0.002,0.002) 

(0.002,0.002) 

(0.002,0.002) 

Wing-Tip  Spacing 

0.025 

0.025 

0.025 

Wing  Index  Dimensions  (ij) 

(41,32) 

(41,32) 

(101,32) 

Store  Index  Dimensions  (ij) 

N/A 

(57,5) 

N/A 

Table  3:  Selected  values  of  CAPTSDv  grid  parameters  (lengths  in  wing  chords,  cw). 


A view  of  grid  G2  in  the  x-y  plane  in  the  neighborhood  of  the  wing/store  configuration  is  shown  in  Figure 
2a  to  illustrate  the  effect  of  grid  clustering  along  the  combined  planform. 

For  the  ENS3DAE  calculations  presented  in  this  paper,  only  the  clean  wing  (to  which  the  store  mass  may 
be  added)  is  considered.  A two-block,  HH-type  grid  is  used  for  this  configuration  and  is  generated  with 
GR1DGEN  VI 3.  The  two  blocks  have  equal  dimensions  and  grid  spacings,  and  correspond  to  the  upper 
and  lower  halves  of  the  computational  space.  It  should  be  noted  that  the  wing-tip  geometry  is  slightly 
modified  from  the  description  given  above  to  close  the  tip  and  simplify  the  grid-generation  process. 
Characteristics  of  the  grid  are  given  in  Table  4. 

Contrary  to  the  convention  shown  in  Figure  la,  in  ENS3DAE,  the  physical  coordinates  are  defined  with  x 
streamwise,  y normal,  and  z spanwise.  However,  consistent  with  CAPTSDv,  the  computational 
coordinates  are  defined  with  i streamwise,  j spanwise,  and  k normal.  Figure  2b  shows  a portion  of  the 
grid  in  the  root  plane.  The  inset  shows  a highly  magnified  view  of  the  leading  edge  and  highlights  a gap 
between  the  upper  and  lower  blocks  equal  to  0.053%  of  the  chord.  A sharp  leading  edge  is  difficult  to 
model  using  an  Euler/Navier-Stokes  solver;  the  gap  helps  by  removing  the  leading  edge  from  the 
computational  domain.  The  leading  edge  is  captured  numerically  through  the  boundary  conditions 
applied  at  points  around  the  physical  leading  edge.  In  the  spanwise  direction,  the  wing  tip  is  modeled  by 
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transitioning  the  geometry  from  the  full  airfoil  thickness  at  the  wing  tip  to  zero  thickness  at  the  first 
spanwise  station  outboard  of  the  wing  tip. 


Parameter 

Value 

Domain  Size  (x,y^) 

(406,400,80) 

Block  Dimensions  (i,j,k) 

(161,71,51) 

Leading-Edge  Spacing 

0.02 

Trailing-Edge  Spacing 

0.02 

Normal  Spacing 

0.0006667 

Wing-Tip  Spacing 

0.05 

Wing  Index  Dimensions 

(81,41) 

Table  4:  Selected  values  of  ENS3DAE  grid  parameters  (lengths  in  feet). 


Modal  Analysis 

Modes  of  the  structural  model  are  computed  with  MSC/NASTRAN  and  then  splined  to  aerodynamic 
surface  grids  (specified  at  z = 0)  with  the  infinite  plate  spline,  as  implemented  by  Smith  et  al.  (1995, 
1996).  The  modes  are  scaled  to  yield  generalized  masses  of  magnitude  1.  Modal  amplitudes  are 
computed  through  time  integration  of  the  generalized  structural  dynamics  equations  to  yield  updated 
approximations  of  wing/store  surface  deflections.  Unless  otherwise  stated,  results  given  in  this  paper  are 
obtained  by  retaining  the  4 modes  of  lowest  frequency  in  the  aeroelastic  analysis  and  by  excluding  in- 
plane modes.  Sets  of  retained  of  modes  are  shown  in  Figure  3 for  two  cases:  clean  wing  with  no  store 
and  wing  with  store  mass  positioned  1.75  feet  forward  of  the  wing  elastic  axis  (see  discussion  below). 
With  the  exception  of  mode  4,  the  modes  in  the  two  sets  are  quite  similar;  slightly  more  torsion  is  evident 
at  the  wing  tip  in  the  modes  associated  with  the  clean  wing. 

Linear  Analysis 

To  determine  a parameter  space  (velocity  and  Mach  number  for  specified  altitude)  where  store-induced 
LCO  possibly  exists,  a linear  flutter  analysis  of  a clear  wing  and  a wing/tip  store  combination  is 
conducted.  The  linear  flutter  speeds  (those  from  linear  aerodynamic  theories)  are  determined  from  data 
calculated  from  the  p-k  method  of  MSC/NASTRAN.  The  flutter  and  divergence  instabilities  can  be 
determined  from  an  inspection  of  calculated  data  in  graphical  form,  the  so-called  V-g  and  K-co  diagrams. 
These  diagrams  are  shown  in  Figures  4a  (clean  wing)  and  4b  (wing  with  tip  store  mass)  for  a selected 
Mach  number  of  0.92.  From  Figure  4a,  the  flutter  speed  of  a clean  wing  is  determined  by  the  first 
crossing  of  one  of  the  modes  from  negative  to  positive  values  of  the  damping  parameter,  g (i.e.,  334 
ft/sec),  and  a corresponding  flutter  frequency  of  2.17  HZ.  Additionally,  shown  in  Figures  4a  is  a 
divergence  instability,  whose  speed  is  determined  by  the  simultaneous  occurrence  of  zero  damping  and 
zero  frequency  for  another  mode  (i.e.,  630  ft/sec). 

The  V-g  diagram  of  the  clean  wing  is  compared  to  the  V-g  diagram  of  the  wing/store  configuration  (mass 
only),  shown  in  Figure  4b,  when  the  store  eg  is  located  1.75  ft  upstream  of  the  wing  ea  (equivalent  to  the 
store  pitch  axis).  By  comparing  the  two  diagrams,  it  is  seen  that  the  flutter  speed  is  increased  to  559  ft/sec 
when  the  store  eg  is  placed  at  this  position  and  that  the  severity  of  the  flutter  instability  of  the  clean  wing 
has  been  reduced  (reflected  by  less  damping).  The  flutter  mode  of  the  wing/store  has  been  converted  into 
a “hump,”  or  lightly  damped  mode.  It  is  speculated  that  the  initiation  of  store-induced  LCO  is  associated, 
in  some  way,  with  hump  modes,  such  that  the  linear  flutter  investigation  defined  a beginning  region  to 
search  for  LCO.  Of  course,  since  this  hump  is  determined  from  linear  aerodynamics,  the  region  of  LCO 
must  be  modified  by  taking  into  consideration  transonic  (nonlinear)  aerodynamics.  This  modification  is 
discussed  in  the  following  sections  for  a determination  of  store-induced  LCO. 


44-8 


Linear  analysis  is  carried  out  for  other  store  mass  positions,  but  not  reported  herein.  These  results  show 
the  reduction  of  the  peak  damping  parameter  with  forward  movement  of  the  store  mass.  The  offset 
position  of  1.75  ft  (upstream  of  the  wing  ea)  is  selected  for  use  in  the  CAPTSDv  calculations  reported 
below,  because  of  the  small  peak  value  of  g attained  with  this  parameter  value. 

Flutter  Boundaries 

Boundaries  of  flutter  and  LCO  onset  are  computed  for  the  clean  wing  configuration  and  for  the  wing  with 
store  mass  (i.e.,  store  not  modeled  aerodynamically).  The  flow  is  assumed  to  be  inviscid.  These  results 
are  compared  to  those  obtained  with  MSC/NASTRAN  using  linear  analysis.  LCO  solutions  are  observed 
at  Mach  numbers  above  Mach  0.9  when  the  store  mass  is  present;  these  cases  will  be  discussed  in  greater 
detail  in  the  next  section.  Flutter  and  LCO  boundaries  are  compared  in  Figure  5a.  For  the  clean  wing, 
CAPTSDv  predicts  a flutter  speed  of  433  ft/sec  at  Mach  0.7,  a value  3.5%  higher  than  that  predicted  by 
MSC/NASTRAN.  With  the  store  mass  included,  CAPTSDv  predicts  a flutter  speed  of  648.5  ft/sec 
(CAPTSDv),  a value  6.8%  higher  than  that  provided  by  MSC/NASTRAN.  At  Mach  0.7,  the 
aerodynamics  are  linear  and  the  reasonable  comparisons  between  CAPTSDv  and  MSC/NASTRAN  are  to 
be  expected.  CAPTSDv  clearly  confirms  that  forward  movement  of  the  store  mass  has  a stabilizing 
effect  on  the  aeroelastic  system  for  Mach  numbers  at  or  below  0.9. 

As  Mach  number  increases  beyond  0.7,  the  clean-wing  flutter  boundary  obtained  with  CAPTSDv 
develops  a transonic  dip  with  a minimum  flutter  speed  (355  ft/sec)  at  about  Mach  0.88.  The  boundary  is 
much  flatter  when  the  store  mass  is  included;  flutter  speed  averages  around  645  ft/sec.  Both  boundaries 
show  a rapid  increase  in  flutter  speed  near  Mach  0.9.  Flowever,  at  selected  Mach  numbers  between  0.90 
and  0.95  (0.91,  0.92,  0.93,  and  0.94)  with  the  store  mass  present,  LCO  solutions  are  observed.  These 
nonlinear  oscillations  are  computed  at  flight  speeds  much  lower  than  the  nominal,  wing/store  flutter 
speed,  and,  for  some  Mach  numbers,  lower  than  the  clean-wing  flutter  speed.  Thus,  the  presence  of  the 
store  destabilizes  the  system  at  higher  Mach  numbers  in  an  adverse  manner,  i.e.,  to  lower  flight  speeds. 

On  the  flutter  boundary  for  the  wing/store  configuration,  two  different  flutter  modes  are  observed.  These 
are  contrasted  in  terms  of  the  computed  lift  and  moment  coefficients  (taken  about  the  leading  edge)  for 
Mach  0.84  and  Mach  0.9,  as  shown  in  Figures  5b  and  5c.  Note  that  unstable  test  points  are  selected 
above  the  flutter  boundary:  U = 750  ft/sec  at  Mach  0.84  (flutter  at  642.5  ft/sec)  and  U = 850  ft/sec  at 
Mach  0.9  (flutter  at  about  825  ft/sec).  Two  different  frequencies  of  divergent  oscillation  are  observed: 
1.90  Hz  (Mach  0.84)  and  9.52  Hz  (Mach  0.90).  These  frequencies  are  somewhat  larger  than  the  natural 
frequencies  corresponding  to  modes  1 (1.69  Hz)  and  3 (9.17  Hz),  respectively.  For  both  flutter  modes, 
the  phase  relationships  between  peak  lift  and  moment  are  the  same  (about  180  degrees  out  of  phase). 
Differences  in  frequency  correlate  well  with  differences  in  modal  participation  between  the  flutter  modes. 
As  shown  in  Figures  5d  and  5e,  response  is  dominated  by  modes  1 and  2 at  Mach  0.84,  whereas  modes  3 
and  4 dominate  the  response  at  Mach  0.9.  The  impact  of  varying  modal  participation  on  wing  shape  is 
described  in  the  next  section. 

A switching  of  flutter  modes  is  perhaps  suggested  by  the  linear  MSC/NASTRAN  results  shown  for  Mach 
0.92  in  Figure  4 (the  results  for  Mach  0.9  are  not  markedly  different).  The  effect  of  nonlinearity  appears 
to  be  large  in  terms  of  stabilizing  the  interaction  between  modes  1 and  2,  predicted  by  MSC/NASTRAN 
to  occur  at  about  560  ft/sec.  At  a larger  frequency,  coalescence  of  modes  3 and  4 is  evident  in  the 
MSC/NASTRAN  results  at  a flight  speed  of  about  800  ft/sec,  a velocity  near  the  flutter  speed  predicted 
by  CAPTSDv. 

Store-Induced  Limit-Cycle  Oscillation 

In  a Mach  number  range  between  0.91  and  0.95,  fully  developed  LCO  states  are  computed  with 
CAPTSDv  for  the  wing/store  configuration,  excluding  store  aerodynamics  (see  Figure  5a).  As  will  be 
shown  later,  the  effect  of  including  store  aerodynamics  is  not  significant,  while  the  effect  of  viscosity  is 
to  reduce  LCO  amplitude.  Two  types  of  LCO  are  observed:  (1)  an  expected  form  involving  significant 
time-periodic  oscillations  of  the  aeroelastic  system  that  will  be  described  first  and  referred  to  as  simply 
LCO,  and  (2)  an  unexpected  form  with  very  small  amplitudes  (~3  orders  of  magnitude  less  in  magnitude) 
that  will  be  described  second  and  referred  to  as  “embryonic”  LCO,  or  ELCO. 
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As  shown  in  Figure  5a,  LCO  is  observed  over  a restricted  range  of  Mach  numbers.  Generally,  LCO 
amplitudes  increase  with  increasing  velocity,  and  for  sufficiently  large  velocities,  computed  oscillations 
become  so  large  that  the  assumptions  of  TSDT  become  invalid.  Attention  is  first  given  to  Mach  0.92, 
where  LCO  is  first  observed  at  U = 390  ft/sec  for  inviscid  flow  (onset  occurring  between  385  and  390 
ft/sec)  and  U = 410  ft/sec  for  viscous  flow  (onset  occurring  between  390  and  410  ft/sec).  It  should  be 
noted  that  the  onset  of  LCO  is  computationally  expensive  to  obtain,  since  very  large  integration  times  are 
necessary  for  the  aeroelastic  system  to  approach  time-asymptotic  behavior  at  values  of  flight  speed  near 
critical.  In  these  calculations,  the  initial  state  is  defined  by  a flow  solution  given  by  a rigid-body 
calculation,  and  the  modal  amplitudes  are  assumed  to  vanish,  except  for  the  first  mode,  which  is  assigned 
an  initial  value  of  0.01. 

The  slow  growth  of  lift  coefficient  to  its  asymptotic  value  is  shown  in  Figure  6a  for  U = 410  ft/sec, 
assuming  inviscid  flow.  LCO  frequency  is  observed  to  be  2.92  Hz.  Time  histories  of  fully  developed  lift 
and  moment  coefficient  are  shown  in  Figure  6b;  the  phase  relationship  between  these  quantities  is  similar 
to  that  observed  for  the  two  flutter  modes  previously  described.  The  computed  frequency  corresponds  to 
the  modal  content  of  the  aeroelastic  response.  As  seen  in  Figure  6c,  mode  2 (a  natural  frequency  of  3.05 
Hz)  dominates  the  response,  with  strong  participation  also  from  mode  1.  At  this  Mach  number  and  flight 
speed,  the  first  two  modes  are  much  more  strongly  coupled  than  in  the  flutter  mode  found  at  low 
transonic  Mach  numbers  (cf.  Mach  0.84),  with  an  associated  increase  in  response  frequency. 

A similar  set  of  plots  of  LCO  response  are  shown  in  Figure  7 for  the  case  of  viscous  flow,  assuming  an 
equivalent  flight  speed  of  410  ft/sec  at  Mach  0.92.  Grid  G3  is  used  for  viscous  calculations,  and  the 
aerodynamics  of  the  store  are  assumed  negligible.  LCO  amplitude  is  observed  to  diminish  through  the 
effects  of  viscosity,  and  frequency  slightly  increases  to  2.99  Hz.  Viscous  simulations  of  LCO  are 
somewhat  stiffer  than  inviscid  computations,  requiring  one  Newton  sub-iterate  per  time  step  for  stable 
calculation. 

Structural  response  associated  with  LCO  is  contrasted  with  that  of  the  two  flutter  modes  through 
visualization  of  wing-tip  motion.  Snapshots  of  instantaneous  chordline  orientation  (lines  connecting  wing 
tip  leading  edge  and  trailing  edge)  are  shown  in  Figure  8 for  each  of  the  three  characteristic  aeroelastic 
responses.  The  low-speed  flutter  mode  is  primarily  first  bending,  which  is  reflected  by  a vertical 
displacement  of  the  wing  tip  at  different  times,  with  introduction  of  only  a slight  incidence  angle  with 
respect  to  the  freestream.  The  flutter  mode  observed  at  Mach  0.9  is  a higher  frequency  mode  involving 
significant  contributions  from  modes  3 and  4.  The  response  of  the  wing  tip  involves  both  pitch  and 
plunge,  with  a motion  akin  to  that  of  a ship  moving  into  ocean  waves  (a  visual  rotation  about  the  mid- 
chord). In  the  situation  of  LCO,  the  response  is  composed  of  first  bending  and  first  torsion  contributions. 
A pitching  motion  dominates  the  resulting  movement  of  the  wing  tip,  with  visual  rotation  about  the  tip 
leading  edge. 

LCO  solutions  of  the  type  described  above  are  not  observed  at  flight  speeds  below  390  ft/sec  for  inviscid 
flow  at  Mach  0.92.  However,  at  speeds  between  340  and  390  ft/sec,  sustained  oscillations  of  very  small 
amplitude  and  less  regular  character  are  computed.  In  this  speed  range,  the  magnitude  of  system 
oscillations  increases  very  slowly  with  increasing  U.  These  ELCO  states  appear  to  be  physical  and  not 
numerical  in  origin,  since  ELCO  formation  is  found  to  be  persistent  to  variation  of  Mach  number  and 
various  numerical  parameters.  However,  ELCO  amplitude  is  sensitive  to  structural  damping  and  the 
numerical  precision  of  the  computation.  When  C,  is  increased  from  the  baseline  value  of  0 to  0.03,  ELCO 
amplitude  is  reduced  by  over  a factor  of  5 at  U = 385  ft/sec,  and  is  found  to  vanish  at  U = 350  ft/sec. 
Also,  increasing  the  baseline  precision  of  the  computation  to  double,  ELCO  amplitude  is  observed  to 
grow  considerably  (but  remaining  at  levels  small  compared  to  LCO). 

A time  history  of  lift  coefficient  is  shown  in  Figure  9 for  U=  385  ft/sec  and  assuming  baseline  values  of 
numerical  parameters:  the  frequency  of  2.94  Hz  is  nearly  identical  to  that  found  during  LCO  at  (7=410 
ft/sec,  while  peak  lift  coefficient  reaches  only  about  5 x 1 O '. 

LCO  Sensitivities 

The  variation  of  LCO  amplitude  with  respect  to  changes  in  flight  speed  at  Mach  0.92  is  computed  for 
three  categories  of  analysis:  inviscid  analysis  of  wing  with  store  mass  (grid  Gl);  inviscid  analysis  of  wing 


44-10 


with  store  modeled  aerodynamically  (grid  G2),  and  viscous  analysis  of  wing  with  store  mass  (grid  G3). 
Results  are  compared  in  Figure  10a  to  show  the  effects  of  varying  the  level  of  modeling  fidelity  within 
the  context  of  transonic  small-disturbance  theory.  Assuming  inviscid  flow,  it  is  observed  that  modelling 
the  store  aerodynamically  has  little  impact  on  the  onset  or  computed  amplitude  of  LCO.  However,  at 
speeds  exceeding  420  ft/sec,  solutions  can  not  be  stably  computed.  In  these  cases  of  increased  wing-tip 
twist,  numerical  destabilization  appears  to  be  a result  of  a very  large,  localized,  pressure  spike  observed 
in  the  region  of  the  juncture  between  the  store  and  the  wing  leading  edge.  This  destabilization  occurs  at 
about  the  speed  for  which  LCO  amplitude  is  predicted  to  grow  rapidly  when  store  aerodynamics  are 
ignored;  peak  lift  coefficient  begins  to  take  a large  jump  at  U = 427  ft/sec  with  this  approximation.  At 
speeds  exceeding  this  value,  the  validity  of  the  computed  inviscid  solutions  is  considered  to  be 
diminished,  owing  to  the  small-disturbance  nature  of  the  methodology.  As  described  above  for  U = 410 
ft/sec,  the  effect  of  viscosity  is  to  diminish  LCO  amplitude.  When  the  boundary  layer  thickness  is 
modelled,  no  large  increase  in  LCO  amplitude  is  observed,  and  peak  lift  coefficient  remains  bounded  by 
0.26  over  the  range  of  speeds  examined,  thus  reducing  the  deformation  of  the  wing  and  extending  the 
speed  range  over  which  the  assumption  of  small  disturbances  is  arguably  satisfied.  Near  the  bifurcation 
point,  at  U = 388  ft/sec,  non-unique  states  are  observed.  One  state,  obtained  with  the  baseline  initial 
conditions  exhibits  very  small-amplitude  ELCO  behavior,  while  the  other  state,  obtained  by  restarting  the 
aeroelastic  solution  from  U = 390  ft/sec,  exhibits  LCO  behavior.  However,  as  velocity  is  reduced  below 
388  ft/sec,  only  ELCO  is  observed,  and  as  velocity  is  increased  above  this  same  speed,  only  LCO  is 
observed.  These  results  are  indicative  of  the  presence  of  a subcritical  bifurcation  just  above  this  speed, 
causing  relatively  large  jumps  in  LCO  amplitude  over  a small  range  of  flight  speeds.  This  phenomenon 
is  explored  further  below  for  Mach  0.93. 

The  variation  of  peak  lift  coefficient  with  velocity  is  found  to  be  more  regular  at  Mach  0.93  than  at  Mach 
0.92.  For  this  higher  Mach  number,  the  effects  of  variation  of  root  angle-of-attack  are  assessed  in  Figure 
10b  using  the  observed  half-differentials  between  maximum  and  minimum  lift  coefficient  (to  separate  the 
dependence  of  oscillation  amplitude  on  velocity  from  the  dependence  of  static  aeroelastic  response  on 
velocity).  When  a0  is  increased  to  2 degrees,  a considerable  decrease  in  LCO  amplitude  is  observed, 
along  with  a significant  delay  in  LCO  onset  to  higher  velocities. 

As  stated  above  for  the  results  computed  at  Mach  0.92  (cf.  Figure  10a),  the  onset  of  LCO  is  subcritical. 
This  is  also  observed  at  Mach  0.93,  as  shown  in  detail  in  Figure  11  (for  the  baseline  condition  of 
vanishing  root  angle-of-attack).  For  this  comparison,  aeroelastic  solutions  are  computed  using  two  kinds 
of  initial  conditions:  the  initial  conditions  described  above,  and  initialisation  of  the  fiowfield  using  a fully 
developed  LCO  solution  obtained  at  a higher  velocity.  It  is  seen  that  the  latter  class  of  initial  condition 
produces  fully  developed  LCO  solutions  for  a small  range  of  velocities  (noted  at  386  ft/sec  and  388 
ft/sec)  at  which  the  former  class  of  initial  condition  does  not  produce  LCO.  While  the  range  of  hysteresis 
is  slight,  the  subcritical  nature  of  the  bifurcation  does  explain  the  rather  large  jumps  in  amplitude 
observed  beyond  the  critical  points. 

Not  surprisingly,  over  the  range  of  Mach  numbers  that  sustain  LCO,  amplitude  of  response  is  found  to  be 
highest  over  the  midsection  of  the  range.  This  is  shown  in  Figure  12  for  a velocity  of  410  ft/sec.  The 
response  near  the  low-Mach  boundary  of  the  LCO  region  is  characteristic  of  the  subcritical  response 
shown  in  Figures  10a  and  11.  However,  at  the  high-Mach  boundary  of  the  LCO  region,  the  computed 
results  are  not  suggestive  of  non-unique  flow  responses. 

Comparisons  with  ENS3DAE 

To  assist  in  the  validation  of  the  transonic  aeroelastic  modeling  employed  in  CAPTSDv,  CAPTSDv 
solutions  (inviscid)  are  compared  to  ENS3DAE  (Euler)  solutions  for  two  clean-wing  calculations: 
prediction  of  static  aeroelastic  response  for  a root  angle-of-attack  of  2 degrees  and  dynamic  response  of  a 
rigid  wing  undergoing  3-Hz  pitch  oscillations  of  'A  degree  in  amplitude.  Both  cases  show  good  agreement 
between  the  Euler  and  TSDT  results,  and  are  reported  here  for  a transonic  condition  (U  = 400  ft/sec  and 
Mach  0.92)  in  terms  of  the  predicted  values  of  lift  and  moment  coefficients.  For  the  case  of  static 
aeroelastic  response,  CL  predictions  agree  within  6%  (0.126  for  ENS3DAE  and  0.134  for  CAPTSDv)  and 
CM  predictions  agree  within  15%  (-0.0687  for  ENS3DAE  and  -0.0809  for  CAPTSDv).  In  the  case  of 
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response  to  forced  pitch  oscillations,  lift  and  moment  coefficients  are  compared  in  Figure  13  as  phase 
portraits,  where  it  is  seen  that  amplitudes  and  phase  relationships  are  in  excellent  agreement. 

Summary  and  Conclusions 

A class  of  limit-cycle  oscillations  was  observed  for  a rectangular  wing  with  tip  store.  These  LCOs 
occurred  at  speeds  lower  than  that  predicted  using  linear  aerodynamics  and  at  speeds  lower  than  that 
computed  for  the  clean-wing  configuration.  The  form  of  the  bifurcation  was  subcritical,  such  that  LCO 
amplitude  jumped  abruptly  as  Mach  number  increased  beyond  a critical  value.  However,  it  was  also 
found  that  as  Mach  number  increased  to  a critical  value  (~0.94-0.95),  LCO  states  could  no  longer  be 
sustained.  For  the  configurations  examined,  the  presence  of  LCO  was  insensitive  to  the  inclusion  of 
store  aerodynamics  in  the  aeroelastic  model.  Also,  the  effect  of  viscosity  was  to  diminish  LCO 
amplitude.  A second  class  of  LCO  solutions  with  small  amplitude  was  observed  that  occurred  over  a 
range  of  speeds  below  critical,  i.e.,  prior  to  the  initiation  of  LCOs  characterized  by  large-amplitude 
aeroelastic  response.  These  states  were  found  to  be  sensitive  to  structural  damping,  such  that  addition  of 
nominal  levels  of  damping  was  sufficient  to  overcome  the  phenomenon. 

The  search  for  LCO  states  was  conducted  in  two  steps.  First,  linear  theory  was  employed  in  the 
identification  of  “hump”  modes,  which  corresponded  to  placement  of  the  store  mass  near  the  wing 
leading  edge.  Such  characteristic  aeroelastic  responses  were  expected  to  point  to  conditions  susceptible 
to  LCO.  Second,  LCO  states  were  computed  using  the  transonic  small-disturbance  theory  algorithm 
CAPTSDv,  assuming  both  inviscid  and  viscous  flow.  This  nonlinear  mathematical  formulation  was 
sufficient  to  capture  properly  weak  shock  formation  and  movement.  The  validation  of  these  LCO  results 
using  higher  fidelity  methods  (i.e.,  Euler  and  Navier-Stokes  via  ENS3DAE)  is  currently  underway,  and 
preliminary  results  are  encouraging.  Also,  static  aeroelastic  and  dynamic  rigid  responses  of  the  clean 
wing  computed  with  CAPTSDv  and  ENS3DAE  are  in  close  agreement.  Additional  work  is  planned  for 
the  study  of  LCO  including  the  presence  of  underwing  stores  and  for  the  study  of  a swept-wing 
configuration.  Further  investigation  is  required  to  yield  a full  explanation  of  LCO  development  for  the 
wing/store  configuration  examined  herein  and  to  better  understand  any  potential  connection  between  the 
observed  nonlinear  LCO  phenomenon  and  the  development  of  hump  modes  in  the  reported  linear 
analysis. 
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Mode  1 
1.97  Hz 


Mode  2 
4 05  Hz 


(a)  Clean  wing 


Mode  3 
9.17  Hz 


(b)  Wing  with  store  mass  (offset  1.75  ft  forward) 


Figure  3:  Splined  modes  retained  in  aeroelastic  analysis  of  Goland  wing  with  and  without  store  mass. 


Damping  Parameter 
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(a)  Clean  wing  damping  and  frequencies 


(b)  Wing/store  damping  and  frequencies 


Figure  4:  Damping  and  frequency  characteristics  predicted  by  linear  aeroelastic  analysis 
(MSC/NASTRAN)  for  the  Goland  wing  with  and  without  store  mass  at  Mach  0.92. 
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Figure  5:  Inviscid  CAPTSDv  aeroelastic  analysis  for  grid  G1  without  store  aerodynamics:  (a)  comparison  of 
flutter  and  LCO  boundaries  computed  with  CAPTSDv  and  MSC/NASTRAN;  (b)  and  (c)  force  and  moment 
coefficient  histories  for  Mach  0.84  ( U=  750  ft/sec)  and  Mach  0.90  (U=  850  ft/sec);  (d)  and  (e)  modal 
amplitude  histories  for  Mach  0.84  (U  = 750  ft/sec)  and  Mach  0.90  (U=  850  ft/sec). 


(b) 


Time,  t 


(C) 


Time,  t 


Figure  6:  Inviscid  CAPTSDv  aeroelastic  analysis  for  grid  G1  without  store  aerodynamics  at  Mach  0.92  and 
(7=41 0 ft/sec:  (a)  time  history  of  lift  coefficient;  (b)  lift  and  moment  coefficients  at  LCO;  (c)  modal 

amplitudes  at  LCO. 
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Figure  7:  Viscous  CAPTSDv  aeroelastic  analysis  for  grid  G3  without  store  aerodynamics  at  Mach  0.92  and 
(7=410  ft/sec:  (a)  time  history  of  lift  coefficient;  (b)  lift  and  moment  coefficients  at  LCO;  (c)  modal 

amplitudes  at  LCO. 
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Figure  10:  Comparison  of  LCO  response  boundaries  for  different  CAPTSDv  analyses:  (a)  inviscid  analysis 
with  (G2)  and  without  (Gl)  inclusion  of  store  aerodynamics,  and  viscous  analysis  (G3)  for  Mach  0.92;  (b) 
inviscid  analysis  without  store  aerodynamics  at  selected  angles  of  attack  for  Mach  0.93. 


375  400  425  450 

Velocity,  U 

Figure  11:  Demonstration  of  subcritical  Hopf  bifurcation  and  non-unique  aeroelastic  behavior  at  Mach  0.93 
(wing/store  model  without  store  aerodynamics).  Two  initial  conditions  are  employed:  baseline  initial 
conditions  (“Perturbed  Rigid  IC”)  and  converged  LCO  solutions  at  higher  velocities  (“Path  Following  1C”). 
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Figure  12:  Variation  of  LCO  response  for  different  Mach  numbers  at  U = 410  ft/sec  (wing/store  model 

without  store  aerodynamics). 


Figure  13:  Comparison  of  CL-CM  phase  portraits  predicted  by  CAPTSDv  and  ENS3DAE  for  forced  pitch 
oscillation  (amplitude  of  Vi  deg  at  3 Hz)  of  the  rigid  wing  at  Mach  0.92  and  U = 400  ft/sec. 


